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Abstract. Spatiotemporal correlations of the two-dimensional spring-block 
(Burridge-Knopoff) model of earthquakes are extensively studied by means 
of numerical computer simulations. The model is found to exhibit either a 
"subcritical" , "supercritical" or "near-critical" behavior, depending on the 
values of the model parameters. Transition between the "subcritical" and "su- 
percritical" regimes is either continuous or discontinuous. Seismic events in 
the "subcritical" regime and those in the "supercritical" regime at larger mag- 
nitudes exhibit universal scaling properties. In the "supercritical" regime, 
eminent spatiotemporal correlations, e.g., a remarkable growth of seismic ac- 
tivity preceding the mainshock, arise in earthquake occurrence, whereas such 
spatiotemporal correlations are significantly suppressed in the "subcritical" 
regime. Seismic activity is generically suppressed just before the mainshock 
in a close vicinity of the epicenter of the upcoming event while it remains 
to be active in the surroundings (the Mogi doughnut). It is also observed that, 
before and after the mainshock, an apparent S-value of the magnitude dis- 
tribution decreases or increases in the "supercritical" or "subcritical" regimes, 
respectively. Such distinct precursory phenomena may open a way to the pre- 
diction of the upcoming large event. 
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1. Introduction 

An earthquake is a stick-slip dynamical instability of a pre-existing fault driven by 
the motion of a tectonic plate [Scholz, 2002;Scholz, 1998]. While an earthquake is a 
complex phenomenon, certain empirical laws such as the Gutenberg-Richter (GR) law 
and the Omori law concerning its statistical properties are known to hold. These are 
both power-laws without any characteristic energy or time scale. Such "critical" features 
of the statistical properties of earthquakes lead to the view that earthquake might be 
a phenomenon of "self-organized criticality (SOC)". This view is opposite to the other 
widely-spread view of "characteristic earthquakes" where an earthquake is regarded to 
possess its characteristic energy or time scale. 

Modeling earthquakes and elucidating its statistical properties have been a fruitful strat- 
egy in earthquake studies. One of the popular model might be the so-called spring-block 
model originally proposed by Burridge and Knopoff (BK) [Burridge and Knopoff, 1967]. 
In this model, an earthquake fault is simulated by an assembly of blocks, each of which is 
connected via the elastic springs to the neighboring blocks and to the moving plate. All 
blocks are subject to the friction force, the source of the nonlinearity in the model, which 
eventually realizes an earthquake-like frictional instability. While the spring-block model 
is obviously a crude model to represent a real earthquake fault, its simplicity enables one 
to study its statistical properties with high precision. 

Carlson, Langer and others [Carlson and Langer, 1989a; Carlson and Langer, 1989b; 
Carlson et al., 1991; Carlson, 1991a; Carlson, 1991b; Carlson et ai, 1994] studied the 
statistical properties of the BK model quite extensively, paying particular attention to 
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the magnitude distribution of earthquake events. Most of these simulations have been 
done for the simplest one-dimensional (ID) version of the model. It was then observed 
that, while smaller events persistently obeyed the GR law, i.e., staying critical or near- 
critical, larger events exhibited a significant deviation from the GR law, being off-critical 
or "characteristic" [Car/son and hanger, 1989a; Carlson and Langer, 1989b; Carlson et 
al, 1991; Carlson, 1991a; Carlson, 1991b; Schmittbuhl et al, 1996]. The spring-block 
model has also been extended in several ways, e.g., taking account of the effect of the 
viscosity [Myers and Langer, 1993; Shaw, 1994; De and Ananthakrisna, 2004], modifying 
the form of the friction force [Myers and Langer, 1993; Shaw, 1995; Cartwright, 1997; De 
and Ananthakrisna, 2004], considering the long-range interactions between blocks [Xia et 
al, 2005, 2007], driving the system only at one end of the system [Vieira, 1992], or by 
incorporating the rate- and state-dependent friction law [Ohmura and Kawamura, 2007]. 

In the previous paper, the present authors studied the statistical properties of the ID 
BK model, focusing on its spatiotemporal correlations [Mori and Kawamura, 2005;2006]. 
This study has revealed several interesting features of the ID BK model. For example, 
preceding the mainshock, the frequency of smaller events is gradually enhanced, whereas, 
just before the mainshock, it is suppressed in a close vicinity of the epicenter of the 
upcoming event (the Mogi doughnut). The time scale of the onset of the doughnut- like 
quiescence depends on the extent of the frictional instability. Furthermore, the apparent 
S-value of the magnitude distribution increases significantly preceding the mainshock 
under certain conditions. 

It should be remembered, however, that real earthquake faults are 2D rather than ID. 
Hence, it is clearly desirable to study the 2D version of the model in order to further clarify 
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the statistical properties of earthquakes. In the present paper, we study spatiotemporal 
correlations of the 2D BK model. The 2D model studied here is intended to represent 
a 2D fault plane itself, where the direction orthogonal to the fault plane is regarded to 
be rigid and not considered explicitly in the model [Carlson, 1991b]. The other possible 
version of the 2D model is the one where the second direction of the model is taken to be 
orthogonal to the fault plane [Myers et al., 1996]. 

Previous calculations on the 2D BK model were performed with main interest in their 
magnitude distribution, while there were very few systematic studies of its spatiotemporal 
correlations. In the present paper, we wish to fill this gap by investigating the spatiotem- 
poral correlations of the 2D BK model, extending our previous study on the ID BK model 
[Mori and Kawamura, 2005; 2006]. Statistical properties of the model are examined here 
for quite a wide range of the model parameters. Concerning the most relevant parame- 
ter representing the extent of the frictional instability a, we study the properties of the 
model with varying the parameter a over all possible range of [0, ooj. Such comprehensive 
studies of the 2D BK model enables us to identify several regions of the parameter space 
characterized by qualitatively different behaviors. 

The present paper is organized as follows. In §2, we introduce the model and ex- 
plain some of the details of our numerical simulation. The results of our simulations are 
presented in §3 and §4. In §3, we show the results of the event-size distribution, and 
construct a "phase diagram" in the parameter space characterizing the model. The prop- 
erties of the mean displacement and the mean number of failed-blocks of each seismic 
event are analyzed. The possible effect of the spatial anisotropy of the elastic param- 
eters is also examined. Then, in section 4, we analyze various types of spatiotemporal 
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correlation functions of seismic events, including the local recurrence-time distribution, 
the time-correlation function of seismic events before and after the mainshock, the time 
development of the spatial correlation function of seismic events before and after the 
mainshock, and the time development of the magnitude distribution function before and 
after the mainshock. Finally, §5 is devoted to summary and discussion. 

2. The model and the method 

Our model is the 2D version of the spring-block BK model, representing a "fault plane" , 
which is taken to be the x — z plane. The plane consists of a 2D square array of blocks 
containing blocks in the x-direction and Nz blocks in the 2;-direction. All Blocks are 
assumed to move only in the x-direction along strike, and are subject to the friction force 
Each block is connected with its four nearest-neighbor blocks via the springs of the 
elastic constant kc, and is also connected to the moving plate via the spring of the elastic 
constant kp. For simplicity, we assume for most part that the spring constants are spatially 
isotropic, i.e., assume that the elastic constant of the spring in the x-direction kcx and 
the one in the 2;-direction kcz are identical, kcx — Kz — K- (Later, the effect of the spatial 
anisotropy of the spring constants kcx 7^ kcz will be studied briefly.) 

The equation of motion for the block at the site (i, j) is given by 

mUij = kpiu't' - Uij) + kc{Ui+i,j + f/jj+i , . 

where m is the mass of a block, t' is the time, Ui^ is the displacement along the x-direction 
of the block (i, j), and v' is the loading rate representing the speed of the plate. 

The equation is made dimensionless in the same way as in [Mori and Kawamura, 2006] , 
i.e., the time t' is measured in units of the characteristic frequency uj — Jkp/m and the 
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displacement Uij in units of the length L = ^{0)/kp, $(0) being a static friction. Then, 
the equation of motion can be written in the dimensionless form as 

+Ui-ij + Uij-i - Auij) - (i){ui), ^ ' 

where t = t'ui is the dimensionless time, Ui^j = Uij/L is the dimensionless displacement 
of the block (i, j), I = ^Jkjkp is the dimensionless stiffness parameter, v — i/' /{Lu) is the 
dimensionless loading rate, and (j){ui) = $(f/i)/$(0) is the dimensionless friction force. 

As the form of the friction force (t){ui) , we use the velocity- weakening friction force which 
is a single-valued function of the velocity. As its explicit functional form, we use the form 
introduced by [Carlson et al, 1991], 

l+2au/(l— (t) ' ' 

where the friction force immediately drops to 1 — u on sliding, and decays toward zero 
with a rate proportional to the parameter a as the velocity increases. The back-slip is 
inhibited by imposing an infinitely large friction for iii < 0, i.e., (f){u < 0) = — oo. 

This friction force represents the velocity-weakening friction force. Real friction force 
is of course more complex, not depending on the velocity alone. Indeed, various types of 
constitutive laws describing rock friction have been proposed in the literature, e.g. the slip- 
weakening friction force [Scholz, 2002;Shaw, 1995;Myers et ai, l996;Scholz, 1998] or the 
rate- and state-dependent friction force [Dietrich, 1979;Ruina, 1983;Scholz, 2002;Scholz, 
1998], etc. Here, in order to investigate the effect of the dimensionality on the spa- 
tiotemporal correlations of the BK model, we assume the simplest version of the velocity- 
weakening friction force. 
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The friction force is characterized by the two parameters, a and a. The former, cr, 
represents an instantaneous drop of the friction force at the onset of the shp, while the 
latter, a, represents the rate of the friction force getting weaker on increasing the sliding 
velocity. The a = case represents the simplest Coulomb friction law where the friction 
force instantaneously drops from the static value 1 to its dynamical value 1 — o" as soon as 
the block begins to slide, and is kept constant on sliding irrespective of the velocity. The 
a — oo case also corresponds to the another Coulomb friction law where the dynamical 
friction immediately drops to zero on sliding. In addition to these frictional-parameters, 
the model possesses one more material parameter, an elastic-parameter I. We try to cover 
a rather wide range of the parameters a, / and a. For example, we take the frictional- 
parameter a in the range a — [0, oo], and systematically examine the a-dependence of 
the results. 

We also assume the loading rate v to be infinitesimally small, and put v — ^ during 
an earthquake event, a very good approximation for real faults \Carlson et al, 1991]. 
Taking this limit ensures that the interval time during successive earthquake events can 
be measured in units of irrespective of particular values of u. 

Numerical details are the same as in [Mori and Kawamura, 2006] . We solve the equation 
of motion (2) by using the Runge-Kutta method of the fourth order, the width of the time 
discretization At being Atu = 10~®. Total number of 10^ ~ 10^ events are generated in 
each run, which are used to perform various averagings. In calculating the observables, 
initial 10^ events are discarded as transients. 

When the model is regarded to represent a strike-slip fault, the z-direction is supposed 
to correspond to the depth direction. Following [Carlson, 1991b], we impose periodic 
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boundary condition in the x-direction and free boundary condition in the 2;-direction. For 
the most part of our calculation, the system size is taken to be — 160 and — 80. 
In some cases, we study several other sizes, — — N, with varying N in the range 
60 < < 480. 

3. Event-size distribution 

In this and following sections, we show the results of our numerical simulations on the 2D 
BK model for various observables. In this section, we analyze the event-size distributions 
of earthquakes. 

3.1. Magnitude distribution 

We define the magnitude of an event, /x, as a logarithm of its moment M, i.e., 



where Aitj^- is the total displacement of the i,j block during a given event and the sum 
is taken over all blocks involved in the event [Carlson et al, 1991]. 

Typical behaviors of the magnitude distribution of the 2D BK model are shown in 
Figs. 1(a) and (b) for the case of I — 3 and a — 0.01, with varying the frictional-parameter 
a. The magnitude distribution i?(/x)d/x represents the rate of events with their magnitudes 
in the range [/i, /i + d/i]. The parameter a is varied in a wide range of < a < oo. 

Fig. 1(a) exhibits R{n) for smaller a. In the range a 5 0.5, i?(/x) bends down rapidly 
at larger magnitudes. Such a behavior of -R(//) is often called "subcritical" . Only small 
events of S 2 occur in this case. For a ^ 0.5, large earthquakes of their magnitudes 
6 S S 8 suddenly appear, while earthquakes of intermediate magnitudes, say, 2 < // < 6, 
remain rather scarce. It means that large and small earthquakes are well separated at 
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a ~ 0.5. Such a sudden appearance of large earthquakes at a = ad — 0.5 coexisting 
with smaller ones has a feature of a "discontinuous transition". Indeed, the maximum 
magnitude of the observed events (within our observation time) Umax suddenly jumps 
from jimax — 4 at q; = 0.4 to Umax — 7 at q; = 0.5, demonstrating the discontinous nature 
of the transition at o; ~ 0.5. In this connection, it should be mentioned that Vasconcelos 
showed that a single-block system exhibited a first-order transition at a = 1 from a stick- 
slip to a creep [Vasconcelos, 1996], whereas this first-order transition becomes apparently 
continuous in the ID many-block system [Vieira et al., 1993; Clancy and Corcoran, 2005]. 
A discontinuous change observed at o; = ccd — 0.5 in the present 2D model may be related 
to the first-order transition of a single-block system, although events observed at a < ad 
in the present 2D model are not really creeps, but rather are stick-slip events of small 
sizes. 

With increasing a further, earthquakes of intermediate magnitudes gradually increase 
their frequency. Fig. 1(b) exhibits R{n) for larger a. In the range of 1 "S a "S 10, 
i?(//) exhibits a pronounced peak structure at a larger magnitude, deviating from the 
GR law at // ^ /i ~ 2, while it exhibits a near straight- line behavior corresponding to 
the GR law at smaller magnitudes fi lS fl- The observed behavior of R{fJ,) agrees with 
the previous findings by Carlson, Langer and collaborators [Carlson and Lan^-er, 1989a; 
Carlson and Langer 1989b; Carlson et al, 1991]. Such a behavior of R{fi) is sometimes 
called "supercritical", since R{ijl) bends up at larger magnitudes (though it eventually 
falls off at still larger magnitudes). According to Carlson et al, smaller events of their 
magnitudes /i ^ jl are localized events, while larger events of their magnitudes /i ^ jl 
are delocalized ones. The existence of a distinct peak structure at a larger magnitude 
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suggests that large earthquakes of ^ /i are more or less characteristic. In order to 
check that the peak structure in R{ii) is not an artifact due to the finite-size effect, the 
system-size dependence of is examined in Fig. 1(c) for the case of a = 3, for N x N 
systems with 60 < < 480. As can be seen from the figure, the characteristic peak at 
a larger magnitude is affected somewhat by the finite-size for N < 240, while the peak 
is almost A^-independent for > 240. This observation indicates that the characteristic 
peak structure of R{iJ,) is an intrinsic property of the bulk system. 

As a increases further, the peak at a larger magnitude becomes less pronounced. At 
a — — 13, R{iJ>) exhibits a near straight-line behavior for a rather wide magnitude 
range, though R{n) falls off rapidly at still larger magnitudes /i ^ 7, indicating that 
the "near-critical" behavior observed for a — ac2 ^13 cannot be regarded as a truly 
asymptotic one. We have checked that this rapid fall-off of R{iJ,) at very large magnitudes 
is a bulk property, not a finite-size effect: See the inset of Fig. 1(c). On further increasing 
a beyond a > — 13, R{iJ>) exhibits again a bending-down "subcritical" behavior. The 
change from the "supercritical" to "subcritical" behavior ai a — ac2 ~ 13 is a continuous 
one, in contrast to the discontinuous change observed at ck = ad — 0.5. 

The behavior of i?(/i) remains qualitatively the same even when the value of / is varied. 
In Figs. 2 and 3, we show R{ii) both for a smaller value of I — 1, and for a larger 
value of /, / = 5, respectively. In either case, the behavior of i?(/i) is qualitatively similar 
to the one shown in Figs.l for the case of Z = 3. In particular, on increasing a, the 
"subcritical" behavior realized at smaller a changes into the "supercritical" behavior via 
a "discontinuous transition" at a = ad, and then, into the "subcritical" behavior realized 
at larger a via a "continuous transition" at a = ac2- 
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In Fig. 4, we summarize the behavior of in the form of a "phase diagram" in the 
frictional-parameter a versus the elastic-parameter I plane for the case of cr = 0.01. The 
phase diagram consists of three distinct regions, two of which are "subcritical" regions and 
one is "supercritical" region. The "phase boundary" between the smaller-a "subcritical" 
region and the "supercritical" region represents a "discontinuous transition", while the 
one between the larger-a "subcritical" region and the "supercritical" region represents 
a "continuous transition". Note that the "phase boundary" is given here by a visual 
inspection of the magnitude distribution and of several other quantities given in 

§3.3 below. No detailed scaling analysis as performed in [Clancy and Corcoran, 2005] has 
been made. The transition between these different "phases", i.e., a "subcritical phase" for 
small a, a "supercritical phase" for intermediate a, and another "subcritical phase" for 
large a, is primarily dictated by the a-value. Since the "continuous transition" in Fig.4 
has a finite slope in the a-l plane, one can also induce the "subcritical" -"supercritical" 
transition by increasing the /-value for a fixed a. A similar transition induced by the 
change in the Z-value was reported also for the ID BK model [Espanol, 1994; Vieira, 1996] 

So far, we have fixed a io a — 0.01. Next, we analyze the cr-dependence of the magnitude 
distribution. In Figs. 5, we show for various values of a both for a larger value of cr, 
(7 = 0.1, and for a smaller value of cr, cr = 0.001 with fixing / = 3. As can be seen from 
the figures, qualitative features of R{ij) remain essentially the same even including its a- 
dependence as in the case of cr = 0.01. In particular, on increasing ct, the system exhibits 
a discontinuous change from the small-a "subcritical" behavior to the intermediate-a 
"supercritical" behavior at ct = ctd, then a continuous change from the "supercritical" 
behavior to the large-o; "subcritical" behavior at o; = ac2- Hence, the a-l phase diagrams 
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for a = 0.1 and a = 0.001 are qualitatively similar to the one of a = 0.01 shown in Fig. 4. 
Concerning the values of ad and ac2, ctd — 0.5 is insensitive to the u-value, while tends 
to decrease with increasing a. Hence, on increasing a, the width of the "supercritical" 
region is gradually narrowed. In fact, we have observed that the "supercritical" region 
vanishes for / = 3 at around a ~ 0.5. Hence, for sufficiently large a close to cr = 1, 
the "supercritical" region is expected to vanish in the a-l phase diagram, only a single 
"subcritical" region being left. 

In Fig. 6, we show -R(Ai) for various values of a with fixing a — 0. Note that this a — 
case corresponds to the case of a constant dynamic friction of strength 1 — a. As can 
be seen from Fig. 5, always exhibits the "subcritical" behavior for all values of a. 

Namely, the constant dynamical friction always gives rise to the "subcritical" behavior. 
In particular, the case (a = 0, cr = 1) corresponding to the vanishing dynamical friction 
is equivalent to the case {a — oo, a). This is fully consistent with our observation that 
the a = oo limit always yields a "subcritical" behavior in Figs. 1-3 and 5. 

3.2. Effect of spatial anisotropy 

So far, we have assumed for simplicity that the elastic-parameter I is spatially isotropic 
within the 2D plane, i.e., have assigned the same Z- values to the spring constants in the 
X- and in the 2;-directions, Ix = h = I- Kumagai et al suggested that, in order to better 
mimic the original continuum crust, the spring constants of the BK model should be set 
spatially anisotropic, i.e., = \/^lz [Kumagai et al., 1999]. Thus, we also analyze in this 
subsection the effect of spatial anisotropy of the elastic-parameter I on 

In Fig. 7, we show for the case of q;=3.5, = \/3, Iz = 1 and a = 0.01, which 

correspond to exactly the same parameter values as used in [Kumagai et al, 1999]. In 
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the same figure, we also show R{n) of the corresponding isotropic system, I = Ix = h = 
(1 + -\/3)/2 (all other parameter are chosen common) for comparison. This value of I is 
chosen to be a mean of the anisotropic /-values given above. Note that the point (a = 3.5, 
I — -\/3)/2, a — 0.01) hes close to the "continuous" phase boundary in Fig.4. As can 
be seen from Fig. 7, i?(/i) of the isotropic system turns out to be quite similar to that of the 
anisotropic system. We thus conclude that the spatial anisotropy of the elastic parameters 
does not cause any qualitative new feature to the properties of the corresponding spatially 
isotropic model. 

3.3. The mean displacement and the mean number of failed-blocks 

The size of an earthquake event is usually measured by its magnitude. Other possible 
measures of event size might be the mean displacement Am, the number of failed-blocks 
N}, (corresponding to the size of rupture zone), and the mean stress-drop. Note that, in 
the BK model, the mean stress-drop of an event is essentially identical with (proportional 
to) the mean displacement of an event. In Figs. 8 and 9, we respectively show the mean 
displacement and the mean number of failed-blocks as a function of the magnitude for 
various values of a. Naturally, all the quantities are increasing functions of the magnitude. 
An interesting observation here is that the data in the "subcritical" regimes are grouped 
into two distinct branches both in Figs. 8(a) and 9(a), each corresponding to the small-a 
and large-a "subcritical" regions of Fig.4. 

As can be seen from Fig. 8(a), the data in the small-a "subcritical" regime (a < ad — 
0.5) lacks events of larger magnitudes and are characterized by smaller displacement, 
while those in the large-o; "subcritical" regime {a > a^^ — 13) are characterized by 
much larger displacement, by about a factor 10^ larger than that in the small- 
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(X "subcritical" regime with the same \x. All the data of the mean displacement Am in 
the "subcritical" regimes collapse, at least approximately, onto these two curves, which 
are both linear in the magnitude with a common slope ~ 0.1. By contrast, the data in 
the "supercritical" regime (ctci < ol < ac'i) exhibit a significantly different behavior. At 
smaller magnitudes 5 5, they exhibit a crossover behavior between these two universal 
"subcritical" curves depending on its a- value: For smaller a close to ctci, the data tend to 
lie closer to the small-o; "subcritical" curve, while for larger a close to aci-, the data tend 
to lie closer to the large-a "subcritical" curve. At larger magnitudes > 5, the data in 
the "supercritical" regime seem to approach the asymptotic straight line again. At larger 
magnitudes, however, finite-size effects are significant. Hence, we show in Fig. 8(b) the 
system-size dependence of the mean displacement at larger magnitudes, ior NxN systems 
with N in the range 60 < < 480. As can be seen from the figure, the "supercritical" 
curve at large magnitudes appear to approach a straight line with a slope ~ 0.1, which, 
however, is not a simple continuation of the large-a "subcritical" curve. 

The existence of the two "subcritical" curves and the asymptotic "supercritical" curve 
at ^ 5 is also clearly visible in Figs. 9 for the magnitude dependence of the mean 
number of failed-blocks N^. The two "subcritical" curves are again both strikingly hnear 
with a common slope ~ 0.9, while the asymptotic "supercritical" curve at // ^ 5 is also 
apparently linear with a slope ~ 0.9. 

In order to get further insights into the nature of the "subcritical" behavior, we show in 
Figs. 10(a) and (b) the magnitude dependence of the mean displacement and of the mean 
number of failed-blocks for various values oia, a — 0.001, 0.01, 0.1, 1, for the case of a = 
which corresponds to a constant dynamical friction of strength 1 — a. As already shown 
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in Fig.6, all these a = cases correspond to the "subcritical" regime. As expected, all the 
data of the mean displacement show linear behaviors with a common slope ~ 0.1, while 
those of the mean number of failed-blocks show linear behaviors with a common slope 
~ 0.9. An interesting observation here is that the data for different values of a are equally 
spaced both in Figs. 10(a) and (b), meaning that, in the case of a constant dynamical 
friction, the mean displacement (the mean number of failed-blocks) is proportional to the 
parameter a (c~^). 

Hence, the 2D BK model exhibits three distinct behaviors depending on the parameters, 
ie., the small-o; "subcritical" behavior at o; < ctd, the "supercritical" behavior at ccd < 
a < ac2 and the large-a "subcritical" behavior at a > ac2- Given its magnitude, the 
small-o; "subcritical" behavior is characterized by smaller displacement, smaller stress- 
drop and larger rupture-zone size (the number of failed-blocks) with more weight on 
smaller magnitudes, while the large-a "subcritical" behavior is characterized by larger 
displacement, larger stress-drop and smaller rupture-zone size with more weight on larger 
magnitudes. Seismic events in the "subcritical" regime exhibit an interesting scaling 
property, i.e., Au oc M°'^ and iV^ oc M^^^, M = exp /i being a seismic moment of an event. 
Seismic events in the "supercritical" regime exhibit, at smaller magnitudes, a crossover 
behavior between the two "subcritical" behaviors, and exhibit the scaling property at 
larger magnitudes. 

4. Spatiotemporal correlations of earthquakes 
4.1. Local recurrence-time distribution 

Earthquake recurrence, i.e., how earthquakes repeat in time, is a question of general 

interest, and is closely related to the issue of criticahty versus periodicity. Characteristic 
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earthquake recurrence would mean the existence of a characteristic time scale for earth- 
quake recurrence, while critical earthquake recurrence would mean the absence of any 
such characteristic time scale. We investigate earthquake recurrence of the 2D BK model 
via the local recurrence-time distribution function. The local recurrence time T is defined 
by the time passed until the next event occurs with its epicenter lying in a vicinity of 
the previous event within distance of r-blocks from the epicenter of the previous event. 
Here, we consider events with their magnitude /i > /ic {l-i'c — 5,3,0), and compute the 
distribution of the local recurrence time T with r — 5. 

In the main panels of Figs. 11 (a) and (b), we show on a log-log plot the computed 
local distribution function P{T) for the case of a = 3 (a), and a = 20 (b), with fixing 
I — 3 and a — 0.01, each corresponding to the "supercritical" and "subcritical" regimes, 
respectively. In the insets, the same data are re-plotted on a semi-logarithmic scale. The 
recurrence time is normalized by its mean T, which is Tu — 14.0, 19.5, 22.5, respectively 
for lie = 0, 3, 5 in the case of Fig.ll(a), and Tu = 0.14, 1.38, 35.5, respectively for Hc — 0, 
3, 5 in the case of Fig. 11(b). 

As can be seen from the main panels of the figure, P{T) exhibits different behaviors 
between in the "supercritical" and "subcritical" regimes: In the "supercritical" regime, as 
shown in Fig. 11(a), P{T) exhibits two peaks (or a peak and a shoulder), one at a shorter 
time T/T ~ 10^^ and the other at a longer time T/T ~ 10^^, whereas in the "subcritical" 
regime, as shown in Fig. 11(b), P{T) with jic — ^ exhibits a single peak only. In either 
case, P{T) exhibits an exponential tail at longer times: See the insets of Figs. 11. Similar 
exponential tail is also observed in real seismicity [Corral, 2004]. 
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The existence of peak structures and an exponential tail in P{T) means the existence 
of characteristic time scales in earthquake recurrence in the 2D BK model. We note that 
similar characteristic behavior of P{T) was also observed in the corresponding ID BK 
model [Mori and Kawamura, 2005; 2006]. 

4.2. Time correlation of events associated with the mainshock 

In real seismicity, large events often accompany foreshocks and aftershocks. In Fig. 12, 
we show the time correlation function between large events (mainshock) and events of 
arbitrary sizes (dominated in number by small events) for various values of the frictional- 
parameter a with fixing / = 3 and a = 0.01. In the figure, we plot the mean number of 
events of arbitrary sizes occurring within 5 blocks from the epicenter of the mainshock 
before {t < 0) and after (t > 0) the mainshock, where the occurrence of the mainshock 
is taken to be the origin of the time t — 0. The average is taken over all large events of 
their magnitudes with /J, > Hc = 5. The number of events are counted here with the time 
bin of Atu = 0.02. 

In the "supercritical" regime, as can be seen from Fig. 12, a remarkable acceleration of 
seismic activity occurs before the mainshock (t < 0). This feature is also observed in the 
corresponding ID BK model [Mori and Kawamura, 2006]. For the case of a = 3 lying 
in the midst of the "supercritical" regime, seismic activity stays low after the mainshock 
{t > 0), and is gradually enhanced as the time passes. This calm period continues during 
about half of the recurrence time. After this calm period, the seismic activity is activated. 
For the cases of a = 1 and 10 near the border of the "subcritical" regime, seismic activity 
after the mainshock still remains relatively high and is gradually suppressed as the time 
passes. We note that such strong time correlations observed in the "supercritical" regime 
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are robust against the change of the magnitude of the mainshock. In particular, even if 
one chooses as a mainshock an event of a smaller magnitude-threshold //c < 5, one still 
gets eminent time correlations quite similar to the ones shown in Fig. 12 as long as one is 
in the "supercritical" regime. 

In the "subcritical" regime, by contrast, the time correlation is almost absent except 
for the suppression of seismicity immediately before the mainshock, as can be seen from 
Fig.l2. 

4.3. Spatial correlations of events before the mainshock 

In this subsection, we examine the time-development of spatial correlations of seismic 
events before the mainshock. In Figs. 13(a) and (b), we show the spatial correlation func- 
tions between the mainshock and the preceding events of arbitrary size, dominated in 
number by small events, for several time periods before the mainshock, for the cases of 
a — 3 {a) and a = 20 (b) with fixing / = 3 and a = 0.01, each corresponding to the 
"supercritical" and "subcritical" regimes, respectively. It represents the conditional prob- 
ability that, provided that a large event with /i > /ic — 5 occurs at a time event of 
arbitrary size occurs at a time t — to < with its epicenter lying within distance r from 
the epicenter of the mainshock. The computed spatial correlation functions are shown as 
a function of r. Insets represent longer-time behaviors. 

In the "supercritical" regime, as shown in Fig. 13(a), the frequency of small events 
are enhanced preceding the mainshock at and around the epicenter of the upcoming 
mainshock. For small enough t, such a cluster of smaller events correlated with the large 
event may be regarded as foreshocks. The eminent spatial correlations observed in the 
"supercritical" regime turns out to be robust against the change of the magnitude of the 
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mainshock. Just before the mainshock, the frequency of smaller events is suppressed in 
a close vicinity of the upcoming mainshock, while it continues to be enhanced in the 
surrounding blocks. The spatial range where the quiescence occurs is narrow, only of a 
few blocks. Such a doughnut-like quiescence phenomenon prior to the mainshock was 
observed also in the ID BK model. In real seismicity, the doughnut-like quiescence has 
been known as the "Mogi doughnut" [Mogi, 1969;Mogi, 1979;5c/iofe, 2002]. 

In the "subcritical" regime, as shown in Fig. 13(b), the seismic acceleration preceding 
the mainshock is hardly discernible, while the doughnut-hke quiescence is still realized. 
The observation is consistent with the result of the time correlations shown in Fig. 12. 

The doughnut-like quiescence is observed both in the "supercritical" and "subcritical" 
regimes in common. The time scale for the appearance of the doughnut-like quiescence 
depends on the a- and /-values. Namely, the time scale of the onset of the doughnut-like 
quiescence tends to be shorter for smaller a or for larger I. The spatial range where the 
quiescence occurs also weakly depends on the parameter a, getting longer for larger a. 
As mentioned, the observed quiescence is quite similar to the one observed in the ID BK 
model. In the case of the ID BK model, the doughnut-like quiescence can be understood 
in terms of characteristic interval-time scales associated with a single-block event [Mori 
and Kawamum, 2006]. In fact, we have confirmed that essentially the same mechanism 
operates even in the present 2D model, i.e., single-block events give rise to the observed 
doughnut-like quiescence phenomena. 

4.4. Spatial correlation of events after the mainshock 

In this subsection, we examine the time-development of spatial correlations of seismic 
events after the mainshock. The calculated spatial correlation functions are shown in 
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Figs.l4(a)-(d), for the cases of the "supercritical" regime, a = 1 (a), a = 3 (b), a = 10 
(c), and in the "subcritical" regimes, a — 20 (d). Other conditions are taken to be the 
same as in Figs. 13. Insets represent shorter-time behaviors. 

As can be seen from these figures, spatial correlations are eminent in the "supercritical" 
regime (Figs. 14 (a)- (c)), but are almost absent in the "subcritical" regime (Fig. 14(d)). 
Details of the spatial correlations observed in the "supercritical" regime, however, depends 
on the parameter values somewhat. As can be seen from Figs. 14(a) and (c), seismic events 
for a — 1 and 10 remain active in the vicinity of the epicenter of the mainshock, which are 
gradually suppressed as the times passes. One may regard this as "aftershock" . However, 
as can be seen from the insets, at shorter time scales, the frequency of events stays almost 
constant in time. Furthermore, even in the time range where seismic activity decays 
monotonically in time, its time-dependence does not obey the Omori law. 

Somewhat different behavior is observed for the case of a = 3 lying in the midst of the 
"supercritical" regime. As can be seen from Fig. 14(b), seismic activity is rather suppressed 
in a vicinity of the epicenter of the mainshock just after the mainshock, staying almost 
constant in time for some period, and is eventually gradually enhanced. Again, the 
eminent spatial correlations observed in the "supercritical" regime is robust against the 
change of the magnitude of the mainshock. 

In any case, aftershocks obeying the Omori-law generically observed in real seismicity 
is not realized in the 2D BK model for any parameter value. This feature is also common 
to the ID BK model [Carlson and Langer, 1989a.; Carlson and Langer, 1989b;Mon and 
Kawamura, 2006]. 
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4.5. Time-dependent magnitude distribution 

In real seismicity, an appreciable change of the S-value of the magnitude distribution 
has been reported preceding large earthquakes: Often a decrease of the -B- value [Suyehiro 
et al, 1964; Jaume and Sykes, 1999; Kawamura, 2006], but sometimes an increase of it 
[Smith, 1981]. Obviously, a possible change in the magnitude distribution preceding the 
mainshock possesses a potential importance in earthquake prediction. 

In Figs. 15, we show the "time-resolved" local magnitude distributions for several time 
periods before the large event for the cases of a = 1 (a), o; = 13 (b) and a — 20 (c) with 
fixing I — 3 and a — 0.01, each corresponding to the "supercritical", "near-critical" and 
"subcritical" regimes, respectively. Only events with their epicenters lying within 5 blocks 
from the upcoming mainshock is counted here. We define the mainshock as a large event 
of A* > A*c = 5. In the "supercritical" case, as can be seen from Fig. 15 (a), an apparent 
S-value describing the smaller magnitude region, S 2, gets smaller as the mainshock 
is approached, i.e., it changes from B ~ 0.89 of the long-time value to ~ 0.65 in the 
time range tu < 0.1 before the mainshock. In the "subcritical" regime, on the other hand, 
an apparent i^-value gets larger as the mainshock is approached, as can be seen from 
Fig. 15(c). In the "near-critical" case, an apparent S-value hardly changes even when the 
mainshock is approached, although the magnitude distribution loses its weight at larger 
magnitudes: See Fig. 15(b). 

In real seismicity, the S-value usually decreases just before the mainshock, although 
the opposite tendency, i.e., an increase of the 5- value, was also reported. Our present 
result might give a hint to understand such a complex change of the S- value observed in 
real earthquakes. The analysis of the JUNEC seismic catalogue, covering earthquakes in 
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Japan area during the period 1985-1998, have revealed that the i?-value is decreased from 
the all-time value B ~ 0.88 to S ~ 0.60 before the mainshock [Kawamura, 2006]. Such a 
change of the -B- value observed in real seismicity turns out to be close even quantitatively 
to the one observed here in Fig. 15 (a). 

In Figs. 16, we also show the time-resolved local magnitude distributions after the main- 
shock for the cases of a = 1 (a), a = 13 (b) and a — 20 (c) with fixing I — 3 and a — 0.01, 
each corresponding to the "supercritical", "near-critical" and "subcritical" regimes, re- 
spectively. Other conditions are the same as those in Figs. 15. Similarly to Figs. 15, an 
apparent S-value after the mainshock decreases or increases appreciably, depending on 
whether one is either in the "supercritical" or "subcritical" regime. Thus, the situation 
before and after the mainshock turn out to be quite similar. The analysis of the JUNEC 
catalogue yields a change of the i?-value after the mainshock from B ~ 0.71 to the all-time 
value B ~ 0.88 [Kawamura, 2006]. Again, this turns out to be rather close to our result of 
Fig. 16(a) where the 5-value has changed from B ~ 0.67 to the all-time value B ~ 0.89. 

5. Summciry and discussion 

Spatiotemporal correlations of the two-dimensional BK model of earthquakes were stud- 
ied in a wide parameter range by means of extensive numerical computer simulations. A 
phase diagram is constructed in the plane of the elastic-parameter I and the frictional- 
parameter a (Fig. 4). The statistical properties of the model turn out to depend on the 
frictional-parameter a most sensitively. The model exhibits a "subcritical" behavior for 
smaller and larger values of a, while it exhibits a "supercritical" behavior for intermediate 
values of a. Transition between the small-a "subcritical" regime and the "supercritical" 
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regime is discontinuous , while the one between the "supercritical" regime and the large-a 
"subcritical" regime is continuous. 

At smaller magnitudes < /i, the magnitude distribution R{ij,) exhibits a behavior 
close to the power law. At larger magnitudes /i^ jl, however, R{iJ,) exhibits a significant 
deviation from the power law. In the "supercritical" regime realized at intermediate values 
of a, R{ii) exhibits a pronounced peak structure at a large magnitude, falling off rapidly at 
still larger magnitudes, while, in the "subcritical" regimes realized at smaller and larger 
values of a, R{ij) exhibits a decay faster than a power law without a peak structure. 
"Near-critical" behavior corresponding to the GR law is observed only near the boundary 
between the "supercritical" and the large-a "subcritical" regimes. In this sense, a true 
power-law distribution is not a generic attribute of the 2D BK model, neither in the ID 
BK model. In other words, the size distribution of the BK model exhibit an off-critical or 
characteristic behavior at larger magnitudes, the GR-like "near-critical" behavior being 
observed only at smaller magnitudes. 

Given its magnitude, the smaller-a "subcritical" behavior is characterized by smaller 
displacement, smaller stress-drop and larger rupture-zone size, while the larger-a "sub- 
critical" behavior is characterized by larger displacement, larger stress-drop and smaller 
rupture-zone size. Seismic events in the "subcritical" regime exhibit an interesting scaling 
property. Namely, the mean displacement at an event Au scales with the seismic moment 
M — expfi as Au ex M°-^, while the mean number of failed-blocks Nb (the rupture-zone 
size) scales as Nf, oc M°-^. Seismic events in the "supercritical" regime exhibit a crossover 
behavior between the two "subcritical" behaviors at smaller magnitudes, but exhibit the 
scaling property at larger magnitudes. 
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We also observed several intriguing precursory phenomena associated with the main- 
shock, particularly in the "supercritical" regime. Preceding the mainshock, the frequency 
of smaller events is gradually enhanced around the epicenter of the upcoming mainshock. 
Such an increase of the seismic activity arises only in the "supercritical" regime, but is 
hardly visible in the "subcritical" regime. Following the seismic acceleration, the fre- 
quency of smaller events is dramatically suppressed just before the mainshock in a close 
vicinity of the epicenter of the upcoming mainshock, while it remains to be active in the 
surroundings (the Mogi doughnut). Such a doughnut- like quiescence is observed both 
in the "subcritical" and "supercritical" regimes in common, and can be accounted for 
as an instability induced by single-block events. On the other hand, the Omori law of 
aftershocks is not realized in the present 2D BK model, neither in the ID BK model. 

An apparent i?-value before the mainshock decreases or increases appreciably, depend- 
ing on whether one is either in the "supercritical" or "subcritical" regime. Almost the 
same change of an apparent 5- value is observed also after the mainshock. Our result might 
give a hint in understanding the observed change of the S-value in real earthquakes. 

The existence of these distinct precursory phenomena in the 2D BK model may open 
a way to the prediction of the time and the position of the upcoming large event. Of 
course, the BK model is a highly simplified model of real earthquake faults, e.g., finite 
discretization of the originally continuum crust being made, highly simplified constitutive 
relation being assumed, etc. Furthermore, although the present 2D model is certainly 
more realistic than the corresponding ID model in that the full 2D fault plane has been 
taken into account, the real earthquake fault might not necessarily be a sharp 2D plane, 
but rather be a more fractal object with complex branch structure: See, e.51., [Kagan, 
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2006]. Nevertheless, we hope that some of the fundamental characteristics of seismic 
events of the 2D BK model as revealed here might provide us with a useful reference in 
understanding the properties of real seismicity. 
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Figure 1. Magnitude distribution R{fi) of earthquake events for the parameters / = 3 
and a = 0.01. Fig. (a) represents R{fi) for smaller values of the frictional-parameter 
< a < 3, while Fig.(b) represents R{fi) for larger values of the frictional-parameter 
(X). The system size is 60 In^Fi^^c^^'the'siie'ciependence of R{fi) is shown 

for the case of a = 3 (main panel) and a = 13 (inset). 
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Figure 2. Magnitude distribution R{fi) of earthquake events for the parameters / = 1 
and cr = 0.01. Fig. (a) represents -R(yu) for smaller values of the frictional-parameter 
< a < 1.5, while Fig.(b) represents R{fi) for larger values of the frictional-parameter 
1.5 < a < oo. The system size is 60 x 60. 
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Figure 3. Magnitude distribution R{fi) of earthquake events for the parameters / = 5 
and a = 0.01. Fig. (a) represents R{fi) for smaller values of the frictional-parameter 
< a < 15, while Fig.(b) represents R{fi) for larger values of the frictional-parameter 
15 < a < oo. The system size is 60 x 60. 
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Figure 4. Phase diagram in the frictional-parameter a versus the elastic-parameter / 



plane. The parameter cr is cr = 0.01. 
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Figure 5. Magnitude distribution R{n) of earthquake events for various values of a in 
the case of the frictional-parameter a = 0.1 (a), and a = 0.001 (b). The elastic parameter 
is / = 3. The system size is 60 x 60. 
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Figure 6. Magnitude distribution R{fi) of earthquake events for various values of 
the frictional-parameter a. The parameter a is fixed to a = 0, which corresponds to a 
constant dynamical friction of strength 1 — a. The system size is 60 x 60. 
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Figure 7. Magnitude distribution R{fi) of earthquake events in the case of spatially 
anisotropic elastic parameters, = \/3 and = 1, compared with the one in the case of 
spatially isotropic elastic parameters with the same mean value, i.e., Ix = h = (1 + V^)/2. 
The parameters a and a are a = 3.5 and a = 0.01. The system size is 160 x 80. 
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Figure 8. (a) The magnitude dependence of tlie mean displacement for various values 
of tlie frictional-parameter a. Tlie parameters / and a are fixed to / = 3 and a = 0.01. 
Small a {a "S 0.5) and large a {a ^ 13) correspond to the "subcritical" regimes, while 
intermediate a (0.5 lS a lS 13) corresponds to the "supercritical" regime. The data in 
the different regimes, i.e., the small-a "subcritical" regime, the "discontinuous transition" 
regime, the "supercritical" regime and the large-a "subcritical" regimes, are indicated by 
different colors in the figures. The system size is 160 x 80. (b) The mean displacement at 
larger magnitudes are shown on an expanded scale for the N x N systems, with varying 
in the range 60 < iV < 480. 
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Figure 9. (a) The magnitude dependence of tlie mean number of failed-blocks of 
eacli event for various values of tlie frictional-parameter a. Tlie parameters I and a are 
fixed to Z = 3 and a = 0.01. Small a (a 5 0.5) and large a {a ^ 13) correspond 
to the "subcritical" regimes, while intermediate a (0.5 S a S 13) corresponds to the 
"supercritical" regime. The data in the different regimes, i.e., the small-a "subcritical" 
regime, the "discontinuous transition" regime, the "supercritical" regime and the large-a 
"subcritical" regimes, are indicated by different colors in the figures. The system size is 
160 X 80. (b) The mean number of failed-blocks at larger magnitudes are shown on an 
expanded scale for the N x N systems, with varying N in the range 60 < < 480. 
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Figure 10. The magnitude dependence of the mean displacement (a), and of the mean 
number of failed-blocks (b) of each event, for various values of the frictional-parameter a 
in the case of a = 0, which corresponds to a constant djTiamical friction of strength 1 — a. 
The parameters I is fixed to / = 3. The system size is 160 x 80. 
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Figure 11. The local recurrence-time distribution of large events with their magnitudes 
greater than /ic (/ic = 5, 3, 0) for a = 3 (a), and for a = 20 (b), each corresponding to the 
"supercritical" and "subcritical" regimes. The parameters / and a are fixed to be / = 3 
and cr = 0.01. The system size is 160 x 80. The recurrence time T is normalized by its 
mean T, which is T = 14.0, 19.5, 22.5, respectively for /i^ = 0, 3, 5 for the case of Fig. (a), 
and T = 0.14, 1.38, 35.5, respectively for /ic = 0, 3, 5 in the case of Fig.(b). The insets 
represent the semi-logarithmic plots including the tail part of the distribution. The tail 
part shows an exponential behavior for both cases of a = 3 and 20. 
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Figure 12. The time correlation function between large events with = 5 (mainshock) 
occurring at time t = and events of arbitrary sizes (dominated in number by small 
events) occurring at time t for various values of a. The parameters I and a are fixed to be 
/ = 3 and cr = 0.01. Events of arbitrary sizes occurring within 5 blocks from the epicenter 
of the mainshock are counted. The negative time t < represents the time before the 
mainshock, while the positive time t > represents the time after the mainshock. The 
average is taken over all large events with its magnitude n > fic = 5. The system size is 
160 X 80. 
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Figure 13. Event frequency preceding tlie large event witli fi > fic = 5 (mainsliock) 
plotted versus r, the distance from the epicenter of the upcoming mainshock for several 
time periods before the mainshock. The parameter a is a = 3 (a), and a = 20 (b) with 
/ = 3 and cr = 0.01, each corresponding to the "supercritical" and "subcritical" regimes. 
The system size is 160 x 80. The insets represent similar plots with longer-time intervals. 
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Figure 14. Event frequency after the large event with /i > /i^ = 5 (mainshock) plotted 

D R A F T ^ ^ , June 18,^ 2008, 6:p2aiii , , r DRAFT 

versus r, the distance from the epicenter of the preceding mainshock, for several time 



periods after the mainshock. The parameter is a = 1 (a), a = 3 (b), a = 10 (c), and 
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Figure 15. Local magnitude distribution before the mainshock with fi > fic = 5, for 
several time periods before the mainshock, for the cases of a = 1 (a), a = 13 (b) and 
a = 20 (c), each corresponding to the "supercritical", "near-critical" and "subcritical" 
regimes, respectively. Events whose epicenter lies within 5 blocks from the epicenter of the 
upcoming mainshock are counted. The parameters / and a are fixed to be Z = 3 and a = 
0.01. The system size is 160 x 80. In Fig.(a)/(c), an apparent -B-value decreases/increases 



before the mainshock, while it stays almost unchanged in Fig.(b). 
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Figure 16. Local magnitude distribution after the mainshock with fi > fic = 5, for 
several time periods after the mainshock, for the cases of a = 1 (a), a = 13 (b) and 
a = 20 (c), each corresponding to the "supercritical", "near-critical" and "subcritical" 
regimes, respectively. Events whose epicenter lies within 5 blocks from the epicenter of the 
preceding mainshock are counted. The parameters / and a are fixed to be / = 3 and a = 
0.01. The system size is 160x80. In Figs.(a)/(c), an apparent 5-value decreases/increases 



after the mainshock, while it stays almost unchanged in Fig.(b). 
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